home *** CD-ROM | disk | FTP | other *** search
/ Monster Media 1996 #15 / Monster Media Number 15 (Monster Media)(July 1996).ISO / prog_c / cuj0696.zip / DWYER.ZIP / LIB / FAC.C < prev    next >
C/C++ Source or Header  |  1995-10-03  |  6KB  |  268 lines
  1. /*                            fac.c
  2.  *
  3.  *    Factorial function
  4.  *
  5.  *
  6.  *
  7.  * SYNOPSIS:
  8.  *
  9.  * double y, fac();
  10.  * int i;
  11.  *
  12.  * y = fac( i );
  13.  *
  14.  *
  15.  *
  16.  * DESCRIPTION:
  17.  *
  18.  * Returns factorial of i  =  1 * 2 * 3 * ... * i.
  19.  * fac(0) = 1.0.
  20.  *
  21.  * Due to machine arithmetic bounds the largest value of
  22.  * i accepted is 33 in DEC arithmetic or 170 in IEEE
  23.  * arithmetic.  Greater values, or negative ones,
  24.  * produce an error message and return MAXNUM.
  25.  *
  26.  *
  27.  *
  28.  * ACCURACY:
  29.  *
  30.  * For i < 34 the values are simply tabulated, and have
  31.  * full machine accuracy.  If i > 55, fac(i) = gamma(i+1);
  32.  * see gamma.c.
  33.  *
  34.  *                      Relative error:
  35.  * arithmetic   domain      peak
  36.  *    IEEE      0, 170    1.4e-15
  37.  *    DEC       0, 33      1.4e-17
  38.  *
  39.  */
  40. /* -------------------------------------------------------------------- */
  41. /*
  42. Cephes Math Library Release 2.0:  April, 1987
  43. Copyright 1984, 1987 by Stephen L. Moshier
  44. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  45. */
  46.  
  47. #include "mconf.h"
  48.  
  49. /* Factorials of integers from 0 through 33 */
  50. #ifdef UNK
  51. static double factbl[] = {
  52.   1.00000000000000000000E0,
  53.   1.00000000000000000000E0,
  54.   2.00000000000000000000E0,
  55.   6.00000000000000000000E0,
  56.   2.40000000000000000000E1,
  57.   1.20000000000000000000E2,
  58.   7.20000000000000000000E2,
  59.   5.04000000000000000000E3,
  60.   4.03200000000000000000E4,
  61.   3.62880000000000000000E5,
  62.   3.62880000000000000000E6,
  63.   3.99168000000000000000E7,
  64.   4.79001600000000000000E8,
  65.   6.22702080000000000000E9,
  66.   8.71782912000000000000E10,
  67.   1.30767436800000000000E12,
  68.   2.09227898880000000000E13,
  69.   3.55687428096000000000E14,
  70.   6.40237370572800000000E15,
  71.   1.21645100408832000000E17,
  72.   2.43290200817664000000E18,
  73.   5.10909421717094400000E19,
  74.   1.12400072777760768000E21,
  75.   2.58520167388849766400E22,
  76.   6.20448401733239439360E23,
  77.   1.55112100433309859840E25,
  78.   4.03291461126605635584E26,
  79.   1.0888869450418352160768E28,
  80.   3.04888344611713860501504E29,
  81.   8.841761993739701954543616E30,
  82.   2.6525285981219105863630848E32,
  83.   8.22283865417792281772556288E33,
  84.   2.6313083693369353016721801216E35,
  85.   8.68331761881188649551819440128E36
  86. };
  87. #define MAXFAC 33
  88. #endif
  89.  
  90. #ifdef DEC
  91. static short factbl[] = {
  92. 0040200,0000000,0000000,0000000,
  93. 0040200,0000000,0000000,0000000,
  94. 0040400,0000000,0000000,0000000,
  95. 0040700,0000000,0000000,0000000,
  96. 0041300,0000000,0000000,0000000,
  97. 0041760,0000000,0000000,0000000,
  98. 0042464,0000000,0000000,0000000,
  99. 0043235,0100000,0000000,0000000,
  100. 0044035,0100000,0000000,0000000,
  101. 0044661,0030000,0000000,0000000,
  102. 0045535,0076000,0000000,0000000,
  103. 0046430,0042500,0000000,0000000,
  104. 0047344,0063740,0000000,0000000,
  105. 0050271,0112146,0000000,0000000,
  106. 0051242,0060731,0040000,0000000,
  107. 0052230,0035673,0126000,0000000,
  108. 0053230,0035673,0126000,0000000,
  109. 0054241,0137567,0063300,0000000,
  110. 0055265,0173546,0051630,0000000,
  111. 0056330,0012711,0101504,0100000,
  112. 0057407,0006635,0171012,0150000,
  113. 0060461,0040737,0046656,0030400,
  114. 0061563,0135223,0005317,0101540,
  115. 0062657,0027031,0127705,0023155,
  116. 0064003,0061223,0041723,0156322,
  117. 0065115,0045006,0014773,0004410,
  118. 0066246,0146044,0172433,0173526,
  119. 0067414,0136077,0027317,0114261,
  120. 0070566,0044556,0110753,0045465,
  121. 0071737,0031214,0032075,0036050,
  122. 0073121,0037543,0070371,0064146,
  123. 0074312,0132550,0052561,0116443,
  124. 0075512,0132550,0052561,0116443,
  125. 0076721,0005423,0114035,0025014
  126. };
  127. #define MAXFAC 33
  128. #endif
  129.  
  130. #ifdef IBMPC
  131. static short factbl[] = {
  132. 0x0000,0x0000,0x0000,0x3ff0,
  133. 0x0000,0x0000,0x0000,0x3ff0,
  134. 0x0000,0x0000,0x0000,0x4000,
  135. 0x0000,0x0000,0x0000,0x4018,
  136. 0x0000,0x0000,0x0000,0x4038,
  137. 0x0000,0x0000,0x0000,0x405e,
  138. 0x0000,0x0000,0x8000,0x4086,
  139. 0x0000,0x0000,0xb000,0x40b3,
  140. 0x0000,0x0000,0xb000,0x40e3,
  141. 0x0000,0x0000,0x2600,0x4116,
  142. 0x0000,0x0000,0xaf80,0x414b,
  143. 0x0000,0x0000,0x08a8,0x4183,
  144. 0x0000,0x0000,0x8cfc,0x41bc,
  145. 0x0000,0xc000,0x328c,0x41f7,
  146. 0x0000,0x2800,0x4c3b,0x4234,
  147. 0x0000,0x7580,0x0777,0x4273,
  148. 0x0000,0x7580,0x0777,0x42b3,
  149. 0x0000,0xecd8,0x37ee,0x42f4,
  150. 0x0000,0xca73,0xbeec,0x4336,
  151. 0x9000,0x3068,0x02b9,0x437b,
  152. 0x5a00,0xbe41,0xe1b3,0x43c0,
  153. 0xc620,0xe9b5,0x283b,0x4406,
  154. 0xf06c,0x6159,0x7752,0x444e,
  155. 0xa4ce,0x35f8,0xe5c3,0x4495,
  156. 0x7b9a,0x687a,0x6c52,0x44e0,
  157. 0x6121,0xc33f,0xa940,0x4529,
  158. 0x7eeb,0x9ea3,0xd984,0x4574,
  159. 0xf316,0xe5d9,0x9787,0x45c1,
  160. 0x6967,0xd23d,0xc92d,0x460e,
  161. 0xa785,0x8687,0xe651,0x465b,
  162. 0x2d0d,0x6e1f,0x27ec,0x46aa,
  163. 0x33a4,0x0aae,0x56ad,0x46f9,
  164. 0x33a4,0x0aae,0x56ad,0x4749,
  165. 0xa541,0x7303,0x2162,0x479a
  166. };
  167. #define MAXFAC 170
  168. #endif
  169.  
  170. #ifdef MIEEE
  171. static short factbl[] = {
  172. 0x3ff0,0x0000,0x0000,0x0000,
  173. 0x3ff0,0x0000,0x0000,0x0000,
  174. 0x4000,0x0000,0x0000,0x0000,
  175. 0x4018,0x0000,0x0000,0x0000,
  176. 0x4038,0x0000,0x0000,0x0000,
  177. 0x405e,0x0000,0x0000,0x0000,
  178. 0x4086,0x8000,0x0000,0x0000,
  179. 0x40b3,0xb000,0x0000,0x0000,
  180. 0x40e3,0xb000,0x0000,0x0000,
  181. 0x4116,0x2600,0x0000,0x0000,
  182. 0x414b,0xaf80,0x0000,0x0000,
  183. 0x4183,0x08a8,0x0000,0x0000,
  184. 0x41bc,0x8cfc,0x0000,0x0000,
  185. 0x41f7,0x328c,0xc000,0x0000,
  186. 0x4234,0x4c3b,0x2800,0x0000,
  187. 0x4273,0x0777,0x7580,0x0000,
  188. 0x42b3,0x0777,0x7580,0x0000,
  189. 0x42f4,0x37ee,0xecd8,0x0000,
  190. 0x4336,0xbeec,0xca73,0x0000,
  191. 0x437b,0x02b9,0x3068,0x9000,
  192. 0x43c0,0xe1b3,0xbe41,0x5a00,
  193. 0x4406,0x283b,0xe9b5,0xc620,
  194. 0x444e,0x7752,0x6159,0xf06c,
  195. 0x4495,0xe5c3,0x35f8,0xa4ce,
  196. 0x44e0,0x6c52,0x687a,0x7b9a,
  197. 0x4529,0xa940,0xc33f,0x6121,
  198. 0x4574,0xd984,0x9ea3,0x7eeb,
  199. 0x45c1,0x9787,0xe5d9,0xf316,
  200. 0x460e,0xc92d,0xd23d,0x6967,
  201. 0x465b,0xe651,0x8687,0xa785,
  202. 0x46aa,0x27ec,0x6e1f,0x2d0d,
  203. 0x46f9,0x56ad,0x0aae,0x33a4,
  204. 0x4749,0x56ad,0x0aae,0x33a4,
  205. 0x479a,0x2162,0x7303,0xa541
  206. };
  207. #define MAXFAC 170
  208. #endif
  209.  
  210.  
  211. extern double MAXNUM;
  212.  
  213. # if defined(__STDC__) || defined(__PROTO__)
  214. double
  215. fac(int i)
  216. # else
  217. double
  218. fac(i)
  219. int    i;
  220. # endif
  221. {
  222. double x, f, n;
  223. int j;
  224. double gamma();
  225.  
  226. if( i < 0 )
  227.     {
  228.     mtherr( "fac", SING );
  229.     return( MAXNUM );
  230.     }
  231.  
  232. if( i > MAXFAC )
  233.     {
  234.     mtherr( "fac", OVERFLOW );
  235.     return( MAXNUM );
  236.     }
  237.  
  238. /* Get answer from table for small i. */
  239. if( i < 34 )
  240.     {
  241. #ifdef UNK
  242.     return( factbl[i] );
  243. #else
  244.     return( *(double *)(&factbl[4*i]) );
  245. #endif
  246.     }
  247. /* Use gamma function for large i. */
  248. if( i > 55 )
  249.     {
  250.     x = i + 1;
  251.     return( gamma(x) );
  252.     }
  253. /* Compute directly for intermediate i. */
  254. n = 34.0;
  255. f = 34.0;
  256. for( j=35; j<=i; j++ )
  257.     {
  258.     n += 1.0;
  259.     f *= n;
  260.     }
  261. #ifdef UNK
  262.     f *= factbl[33];
  263. #else
  264.     f *= *(double *)(&factbl[4*33]);
  265. #endif
  266. return( f );
  267. }
  268.